
# -----------------------------------------------------------
rm(list=ls())
library(plm)
library(gmm)
library(spatialreg)
library(spdep)
library(splm)
library(AER)

# -----------------------------------------------------------
# Define a function to calculate the spillover productivity
swfn <- function(w) {
  sw <- c()
  for(j in 1:N) {
    sw <- c(sw, sum(w[-j])/(N-1))
  }
  return(sw)
}
# -----------------------------------------------------------
# The following function re-arrange a matrix into a vector
long <- function(X) {
  X2 <- c()
  for(i in 1:nrow(X)) {
    X2 <- c(X2, X[i,])
  }
  return(X2)
}
# -----------------------------------------------------------

# ========================================================
# DGP

set.seed(123)

# Parameters
N <- 100; T <- 10; R <- 1000
d0 <- 0.2; d1 <- 0.55; d2 <- 0; d3 <- 0    # productivity parameters
ak <- 0.25; am <- 0.65    # production function parameters
b1 <- 0.8; b2 <- 0.1    # I policy function parameters
r0 <- 0.01; r1 <- 0.6; r2 <- 0.3     # G parameters
de <- c(0.05, 0.075, 0.1, 0.125, 0.15)   # depreciation rate

# Empty objects
INITIAL <- DFDK <- DFDM <- c()   # matrixes for the estimated elasticities
O1a <- O1b <- O2a <- O2b <- O3a <- O3b <- O4a <- O4b <- O5a <- O5b <- O6a <- O6b <- O7a <- O7b <- O8a <- O8b <- c()
G1a <- G1b <- G2a <- G2b <- G3a <- G3b <- G4a <- G4b <- G5a <- G5b <- G6a <- G6b <- G7a <- G7b <- G8a <- G8b <- c()
SO1a <- SO1b <- SO2a <- SO2b <- SO3a <- SO3b <- SO4a <- SO4b <- c()
SG3a <- SG1a <- SG4a <- SG2a <- c()

O1ap <- O1bp <- O2ap <- O2bp <- O3ap <- O3bp <- O4ap <- O4bp <- O5ap <- O5bp <- O6ap <- O6bp <- O7ap <- O7bp <- O8ap <- O8bp <- c()
G1ap <- G1bp <- G2ap <- G2bp <- G3ap <- G3bp <- G4ap <- G4bp <- G5ap <- G5bp <- G6ap <- G6bp <- G7ap <- G7bp <- G8ap <- G8bp <- c()
SO1ap <- SO1bp <- SO2ap <- SO2bp <- SO3ap <- SO3bp <- SO4ap <- SO4bp <- c()
SG3ap <- SG1ap <- SG4ap <- SG2ap <- c()

# spatial weight matrix
NB <- matrix(1/(N-1), N, N); diag(NB) <- 0; 

# Start the loop
for(ii in 1:R) {
  
  # initial values
  w0 <- runif(N, 1, 3)
  K0 <- runif(N, 10, 200)
  G0 <- runif(N, 0, 1); 
  I0 <- K0^b1*exp(b2*w0)
  M0 <- (am*K0^ak*exp(w0))^(1/(1-am))

  # errors
  eta <- matrix(rnorm(N*T, 0, sqrt(0.07)), N, T)
  xi <- matrix(rnorm(N*T, 0, sqrt(0.04)), N, T)   
  eps <- matrix(rnorm(N*T, 0, sqrt(0.01)), N, T)
  
  # Evolution processes
  # We will first generate state variables, and then control variables whose values
  # depend on state variables.
  Y <- K <- M <- I <- W <- G <- matrix(NA, N, T)
  
  for(i in 1:T) {
    if(i==1) {
      W[,i] <- d0 + d1*w0 + d2*swfn(w0) +d3*G0 + xi[,i]
      K[,i] <- (1-de)*K0 + I0
      G[,i] <- r0 + r1*G0 + r2*W[,i] + eps[,i]; 
      I[,i] <- K[,i]^b1*exp(b2*W[,i])
      M[,i] <- (am*K[,i]^ak*exp(W[,i]))^(1/(1-am))
      Y[,i] <- K[,i]^ak*M[,i]^am*exp(W[,i]+eta[,i])
    }
    else {
      W[,i] <- d0 + d1*W[,(i-1)] + d2*swfn(W[,(i-1)]) +d3*G[,(i-1)] + xi[,i]
      K[,i] <- (1-de)*K[,(i-1)] + I[,(i-1)]
      G[,i] <- r0 + r1*G[,(i-1)] + r2*W[,i] + eps[,i]; 
      I[,i] <- K[,i]^b1*exp(b2*W[,i])
      M[,i] <- (am*K[,i]^ak*exp(W[,i]))^(1/(1-am))
      Y[,i] <- K[,i]^ak*M[,i]^am*exp(W[,i]+eta[,i])
    }
  }
  
  pm <- mean(exp(eta))
  data <- data.frame(id=rep(1:N, each = T), year=rep(1:T, N), Y=long(Y), K=long(K), M=long(M), G=long(G), sm=pm*long(M)/long(Y) )
  
  # ----------------------------------------------
  # Start the estimation
  # Step 1
  h.lnbe <- mean(log(data$sm)); h.eta <- - log(data$sm) + h.lnbe; h.e <- mean(exp(h.eta))
  h.betam <- exp(h.lnbe)/h.e
  data$h.eta <- h.eta
  
  # ----------------------------------------------
  data <- pdata.frame(data)
  data$K1 <- lag(data$K); data$M1 <- lag(data$M); data$G1 <- lag(data$G); data$G2 <- lag(data$G, 2)

  data$y <- log(data$Y); data$k <- log(data$K); data$m <- log(data$M)
  data$k1 <- log(data$K1); data$m1 <- log(data$M1);
  data <- data[!is.na(data$k1),]
  data <- as.data.frame(data)
  #  print(cor(data$k, data$k1)); print(cor(data$m, data$m1)); 
  print(ii)
  
  # ----------------------------------------------
  # Step 2
  # Define function "phi" and "h"
  phi <- function(beta) {
    (1-h.betam)*data$m1 - h.lnbe - log(1/pm) - beta[1]*data$k1 
  }
  
  # Optimization
  RHO.temp <- OBJ <- BK <- c()
  for(ik in seq(0.01, 1, by=0.01)) {
    y.star <- data$y - h.betam*data$m - ik*data$k 
    p1 <- as.vector(phi(ik))
    var.poly <- poly(cbind(p1), degree = 2, raw = TRUE)
    mod1 <- lm(y.star~var.poly)
    RHO.temp <- cbind(RHO.temp, mod1$coeff)
    OBJ <- cbind(OBJ, sum(mod1$residuals^2))
    BK <- cbind(BK, ik)
  }
  target <- which.min(OBJ)
  
  # results - elas
  h.betak <- BK[target]
  DFDK <- cbind(DFDK, h.betak); DFDM <- cbind(DFDM, h.betam)
  
  # Post-reg
  # Generate variables
  data$w <- data$y - h.betak*data$k - h.betam*data$m - data$h.eta
  w1 <- as.vector(phi(c(h.betak)))
  
  for(i in 1:nrow(data)) {
    s1.t <- (as.numeric(data$year)==as.numeric(data$year[i])); s1.t[i] <- 0
    s1.b <- sum(s1.t)
    s1 <- if(s1.b>0) s1.t/s1.b else rep(0, length(s1.t));
    data$sG[i] <- sum(s1*data$G); data$sG1[i] <- sum(s1*data$G1); data$sG2[i] <- sum(s1*data$G2,na.rm = TRUE)
#    data$sL[i] <- sum(s1*(data$G>0)); data$sL1[i] <- sum(s1*(data$G1>0)); data$sL2[i] <- sum(s1*(data$G2>0), na.rm = TRUE)
    data$sw[i] <- sum(s1*data$w); data$sw1[i] <- sum(s1*w1)
  }
  
  pval <- function(x) {summary(x)[["coefficients"]][,"Pr(>|t|)"]} 
    
  ols <- ivreg(w~sG|sG1, data = data); O1a <- cbind(O1a, ols$coefficients["sG"]); O1ap <- cbind(O1ap, pval(ols)["sG"])
  ols <- ivreg(w~G+sG|G1+sG1, data = data); O2a <- cbind(O2a, ols$coefficients[c("G","sG")]); O2ap <- cbind(O2ap, pval(ols)[c("G","sG")])
  ols <- lm(w~sG1, data = data); O3a <- cbind(O3a, ols$coefficients["sG1"]); O3ap <- cbind(O3ap, pval(ols)["sG1"])
  ols <- lm(w~G1+sG1, data = data); O4a <- cbind(O4a, ols$coefficients[c("G1","sG1")]); O4ap <- cbind(O4ap, pval(ols)[c("G1","sG1")])
  ols <- lm(w~sG2, data = data); O5a <- cbind(O5a, ols$coefficients["sG2"]); O5ap <- cbind(O5ap, pval(ols)["sG2"])
  ols <- lm(w~G2+sG2, data = data); O6a <- cbind(O6a, ols$coefficients[c("G2","sG2")]); O6ap <- cbind(O6ap, pval(ols)[c("G2","sG2")])
  ols <- lm(w~G1+sG2, data = data); O7a <- cbind(O7a, ols$coefficients[c("G1","sG2")]); O7ap <- cbind(O7ap, pval(ols)[c("G1","sG2")])
  rm(ols)
  
  pval <- function(x) {summary(x)[["coefficients"]][,"Pr(>|t|)"]} 
  
  sols <- ivreg(w~sw|sw1, data = data); SO1a <- cbind(SO1a, sols$coefficients["sw"]); SO1ap <- cbind(SO1ap, pval(sols)["sw"])
  sols <- ivreg(w~G+sw|G1+sw1, data = data); SO2a <- cbind(SO2a, c(sols$coefficients["G"], sols$coefficients["sw"])); SO2ap <- cbind(SO2ap, c(pval(sols)["G"], pval(sols)["sw"]))
  sols <- lm(w~sw1, data = data); SO3a <- cbind(SO3a, sols$coefficients["sw1"]); SO3ap <- cbind(SO3ap, pval(sols)["sw1"])
  sols <- lm(w~G1+sw1, data = data); SO4a <- cbind(SO4a, sols$coefficients[c("G1","sw1")]); SO4ap <- cbind(SO4ap, pval(sols)[c("G1","sw1")])
  rm(sols)
  
}

# k coefficient
rmse_k <- function(x, truevalue=0.25) sqrt(mean((x-truevalue)^2))
mad_k <- function(x, truevalue=0.25) mean(abs(x-truevalue))
n11 <- mean(DFDK); n12 <- median(DFDK); n13 <- sd(DFDK); n14 <- rmse_k(DFDK); n15 <- mad_k(DFDK)
c(n11, n12, n13, n14, n15)

# -----------------------------------------
# regression stat - mean with s.d.

se <- function(x) { paste0("(", format(round(sd(x), 5), nsmall = 5), ")") }
mean2 <- function(x) {format(round(mean(x), 5), nsmall = 5)}

Table <- matrix(NA, 16, 14)
Table[3,1] <- mean2(O1a); Table[4,1] <- se(O1a)
Table[1,3] <- mean2(O2a[1,]); Table[2,3] <- se(O2a[1,]); Table[3,3] <- mean2(O2a[2,]); Table[4,3] <- se(O2a[2,]); 
Table[9,5] <- mean2(O3a); Table[10,5] <- se(O3a) 
Table[7,7] <- mean2(O4a[1,]); Table[8,7] <- se(O4a[1,]); Table[9,7] <- mean2(O4a[2,]); Table[10,7] <- se(O4a[2,]); 
Table[15,9] <- mean2(O5a); Table[16,9] <- se(O5a) 
Table[13,11] <- mean2(O6a[1,]); Table[14,11] <- se(O6a[1,]); Table[15,11] <- mean2(O6a[2,]); Table[16,11] <- se(O6a[2,]); 
Table[7,13] <- mean2(O7a[1,]); Table[8,13] <- se(O7a[1,]); Table[15,13] <- mean2(O7a[2,]); Table[16,13] <- se(O7a[2,]); 
Table <- rbind(rep(1:7, each=2), rep(c("Coeff",""), times=7), Table)
print(Table, quote=FALSE, na.print = "")

Table <- matrix(NA, 10, 8)
Table[3,1] <- mean2(SO1a); Table[4,1] <- se(SO1a)
Table[1,3] <- mean2(SO2a[1,]); Table[2,3] <- se(SO2a[1,]); Table[3,3] <- mean2(SO2a[2,]); Table[4,3] <- se(SO2a[2,])
Table[9,5] <- mean2(SO3a); Table[10,5] <- se(SO3a)
Table[7,7] <- mean2(SO4a[1,]); Table[8,7] <- se(SO4a[1,]); Table[9,7] <- mean2(SO4a[2,]); Table[10,7] <- se(SO4a[2,])
Table <- rbind(rep(1:4, each=2), rep(c("Coeff"," "), times=4), Table)
print(Table, quote=FALSE, na.print = "")

# -----------------------------------------
# regression stat - significance from 0

sig.0 <- function(x,y) format(round(mean((y>0.05)), 5), nsmall = 5)   

Table <- matrix(NA, 16, 14)
Table[3,1] <- sig.0(O1a,O1ap); 
Table[1,3] <- sig.0(O2a[1,],O2ap[1,]);  Table[3,3] <- sig.0(O2a[2,],O2ap[2,]); 
Table[9,5] <- sig.0(O3a,O3ap); 
Table[7,7] <- sig.0(O4a[1,],O4ap[1,]);  Table[9,7] <- sig.0(O4a[2,],O4ap[2,]);
Table[15,9] <- sig.0(O5a,O5ap); 
Table[13,11] <- sig.0(O6a[1,],O6ap[1,]);  Table[15,11] <- sig.0(O6a[2,],O6ap[2,]); 
Table[7,13] <- sig.0(O7a[1,],O7ap[1,]);  Table[15,13] <- sig.0(O7a[2,],O7ap[2,]);  
Table <- rbind(rep(1:7, each=2), rep(c("1-prob",""), times=7), Table)
print(Table, quote=FALSE, na.print = "")

Table <- matrix(NA, 10, 8)
Table[3,1] <- sig.0(SO1a,SO1ap); 
Table[1,3] <- sig.0(SO2a[1,],SO2ap[1,]);  Table[3,3] <- sig.0(SO2a[2,],SO2ap[2,]); 
Table[9,5] <- sig.0(SO3a,SO3ap); 
Table[7,7] <- sig.0(SO4a[1,],SO4ap[1,]);  Table[9,7] <- sig.0(SO4a[2,],SO4ap[2,]); 
Table <- rbind(rep(1:4, each=2), rep(c("1-prob"," "), times=4), Table)
print(Table, quote=FALSE, na.print = "")

# -----------------------------------------
# regression stat - RMSE

rmse <- function(x, truevalue=0) sqrt(mean((x-truevalue)^2))
mad <- function(x, truevalue=0) mean(abs(x-truevalue))

Table <- matrix(NA, 16, 14)
Table[3,1] <- rmse(O1a)
Table[1,3] <- rmse(O2a[1,]);  Table[3,3] <- rmse(O2a[2,]); 
Table[9,5] <- rmse(O3a); 
Table[7,7] <- rmse(O4a[1,]);  Table[9,7] <- rmse(O4a[2,]); 
Table[15,9] <- rmse(O5a); 
Table[13,11] <- rmse(O6a[1,]);  Table[15,11] <- rmse(O6a[2,]); 
Table[7,13] <- rmse(O7a[1,]);  Table[15,13] <- rmse(O7a[2,]); 
Table <- rbind(rep(1:7, each=2), rep(c("RMSE"," "), times=7), Table)
print(Table, quote=FALSE, na.print = "")


Table <- matrix(NA, 10, 8)
Table[3,1] <- rmse(SO1a); 
Table[1,3] <- rmse(SO2a[1,]);Table[3,3] <- rmse(SO2a[2,]); 
Table[9,5] <- rmse(SO3a); 
Table[7,7] <- rmse(SO4a[1,]); Table[9,7] <- rmse(SO4a[2,]); 
Table <- rbind(rep(1:4, each=2), rep(c("RMSE"," "), times=4), Table)
print(Table, quote=FALSE, na.print = "")


